Method of cancelling ghosts from NMR images

ABSTRACT

A method of cancelling ghosts from NMR images. The method involves estimating a phase difference function Δ (n 1 , n 2 ) and using that function to solve a linear system of equations to find the magnitudes of the true object densities at the true image and ghost locations x(n 2 ,n 2 ) and x(n 1 ,n 2  +N/2), respectively, where the image has dimensions N×N s . Experimental values of Δ (n 1 , n 2 ) for a variety of objects indicate that its variation along n 1  is considerably larger than along n 2 . Thus, for each column n 1 , the phase difference function Δ (n 1 , n 2 ) can be modelled as a one-dimensional function of n 2  with two parameters α (n 1 ) and β (n 1 ), which are estimated from the pixels in the 2-D FFT processed reconstructed image Y(n 1 ,n 2 ). These parameters are then used to estimate Δ (n 1 , n 2 ), which is ultimately used to de-ghost the image.

BACKGROUND OF THE INVENTION

The present invention relates to a method for improving reconstructedNMR images and, more specifically, to a method for cancelling ghostsfrom NMR images.

Under ideal conditions, reconstructed NMR images must be positive andreal. This is because, in NMR experiments, the observed signal isFourier transform of density distribution of the object underconsideration, which by definition is a real positive quantity. Inpractice however, even the most straightforward ways of scanning thek_(x) -k_(y) space, (e.g., row by row or column by column scanningwithout date reversal), result in complex images. This is partly due toa shift in the data from the notional origin in k space. If readoutgradient is constant and uniform time samples of data are inverseFourier transformed to generate the image, delay in the time datatranslates into a linear phase shift in the reconstructed image. Thesephase shifts can be easily determined and eliminated, since they do notaffect the magnitude of the reconstructed images.

When NMR data is obtained by scanning the k_(x) -k_(y) space, with datareversal on alternate y lines (where y is the horizontal axis acrosswhich readout occurs), time delays between the start of data acquisitionand the start of the readout pulse are different for even and odd lines.The effect of this on the image manifests itself as a ghost separated byhalf the image size. Under these conditions, if readout gradient isconstant and data is sampled uniformly in time, then the ghost image canbe entirely removed by a first-order phase difference between odd andeven lines. However, when the readout gradient is sinusoidal and theimage is reconstructed by inverse Fourier transforming non-uniformsamples (see, e.g., the method disclosed in U.S. Pat. No. 4,740,748) thedifference between even and odd line delays degrades introduces ghostsand thus the quality of the resolution.

This, together with asymmetry of the sinusoidal readout gradient foreven and odd lines can be modeled by multiplying even and odd parts ofthe NMR image by two separate phase functions φ (n₁,n₂) and θ (n₁,n₂).More specifically, if x(n₁,n₂) denotes the true density distribution ofthe object under consideration, and Y(n₁,n₂) denotes the 2-D inversediscrete Fourier transform of the time data (i.e., the reconstructedghosted image), then the even and odd parts of the observed image,Y_(even) (n₁,n₂) and Y_(odd) (n₁,n₂)can be modeled as: ##EQU1## wherethe dimensions of the reconstructed image are N×N_(s), and the number ofechoes is N. (Note that if the sinusoidal readout y gradient wasidentical for even and odd lines, then the even and odd phase functionswould have been identical. That is, φ (n₁,n₂).tbd.θ (n₁,n₂)).

Having modeled the ghosted image, the objective can be stated asestimation of the true object density distribution x(n₁,n₂) from theobserved ghosted image Y(n₁,n₂).

From equations (1) and (2), it is clear that if the phase functions φ(n₁,n₂) and θ(n₁,n₂) are known for all values of n₁ and n₂, thenx(n₁,n₂) and x(n₁,n₂ +N/2) can be determined from Y_(even) (n₁,n₂)Y_(odd) (n₁,n₂) by solving a 2×2 linear system of equations. Inpractice, the difference between φ (n₁,n₂) and θ (n₁,n₂) can bedetermined experimentally by placing a test object in the upper andlower half of the field of view (FOV) and measuring the differencebetween even and odd parts of the resulting images. More specifically,when the object is in the upper half of the FOV, by definition:

    x(n.sub.1,n.sub.2 +N/2)=0

Substituting this into equation (1) and (2): ##EQU2##

The phase difference between Y_(even) (n₁,n₂) Y_(odd) (n₁,n₂) can beused to obtain ##EQU3##

Similarly, by placing the object in the lower half of the FOV: ##EQU4##

The phase difference between Y_(even) (n₁,n₂) and Y_(odd) (n₁,n₂) can beused to obtain ##EQU5##

Thus, experimental values of Δ (n₁,n₂) and Δ (n₁,n₂ +N/2) can be used todetermine A(n₁,n₂) and B(n₁,n₂) by solving the above linear system ofequations. Once A and B are determined, their magnitudes can be used tofind x(n₁,n₂) and (n₁,n₂ +N/2) respectively.

Experimentally, there are two major drawbacks with the above approach.The first drawback has to do with the fact that the phase differencefunction is a function of the parameters for the NMR experiments. Someof these parameters are the strength of the x, y and z gradients, andthe static magnetic field or the RF. Therefore, to be able to apply thismethod successfully, a different look up table is needed for differentexperimental set ups. The second drawback has to do with the fact thatthe phase difference function Δ (n₁,n₂) is somewhat object dependent.More specifically, although the general shape of Δ (n₁,n₂) does not varydrastically from one object to the next, the change is large enough tointroduce considerable amount of ghosts. The third drawback of the aboveapproach has to do with the fact that obtaining the phase differencefunction Δ (n₁,n₂) of a test object for all values of n₁ and n₂ is anon-trivial task from an experimental point of view. This has to do withfactors such as susceptibility effects.

In short, it has been found that the performance of the above scheme isinadequate for most ghosted images. Accordingly, it would be highlydesirable to process NMR signals using a method of ghost correction inthe form of an algorithm which is automatic and does not require a lookup table.

SUMMARY OF THE INVENTION

The object of the present invention is therefore to provide a method forautomatically eliminating ghosts in NMR signals resulting from thedifference between even and odd line delays in a traversal of k-spaceusing a sinusoidal readout gradient without using a look-up table.

The present invention achieves the foregoing objective by a methodinvolving the steps of:

(a) taking a two dimensional inverse Fourier transform of the time datato obtain the ghosted image Y(n₁,n₂);

(b) computing the signal energy for each column using ##EQU6##

(c) discarding the columns whose signal energy level are below apredetermined threshold;

(d) estimating α (n₁) and β (n₁) for each remaining column of data; i.e.n₁ =0, . . . N_(s) -1, by:

(i) finding the phase difference function Δ (n₁,n₂) for all ghostingpixels of the column; and

(ii) solving the following simultaneous equations to find linear leastsquare estimates of α (n₁) and β (n₁): ##EQU7##

(e) using α (n₁) and β (n₁) in the above equation to find the phasedifference Δ (n₁,n₂) for 0≦n₂ <N; and

(f) using Δ (n₁,n₂) in equation (11) to find A(n₁,n₂) and B(n₁,n₂) for0≦n₂ <N.

From equations (9) and (10), the true object density distribution at(n₁,n₂) and (n₁,n₂ +N/2) are found by taking magnitude of A(n₁,n₂) andB(n₁,n₂), respectively.

BRIEF DESCRIPTION OF THE DRAWINGS

These and other features and advantages of the present invention willbecome apparent when the following text is read in conjunction with theaccompanying drawings, in which:

FIG. 1 shows a flow diagram of the ghost correction method of thepresent invention;

FIG. 2 depicts a reconstructed image with ghosts.

FIG. 3 graphically depicts the phase difference function; and

FIGS. 4a-11b show before/after examples of NMR images processed with thealgorithm of the present invention.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT

Referring first to FIG. 1, a simplified block diagram of the method ofthe invention is shown. First, the raw NMR signal is converted to animage by performing a two dimensional inverse Fast Fourier Transform(FFT) 2. The resultant image has ghosts, which are eliminated byprocessing the image with the ghost elimination algorithm 4 of thepresent invention.

FIG. 2 depicts an image with ghosts. As shown therein, a bright line 6represents the true image, while a less bright line 8 is a ghost of trueimage 6. Assuming that the traversal through k-space involves thesampling of 128 lines, the ghost will be separated by the true image by64 lines, i.e. half the image size.

FIG. 3 graphically depicts the phase difference function Δ (n₁,n₂). Asshown therein, the phase difference function is appropriatelysymmetrical along the center of the n₂ axis. The phase differencefunction is defined by the offset α and slope β of a straight line 10obtained by taking the least square estimation of a number of points;certain points 12,14 are eliminated (ignored) as being out-of-rangesince Δ (n₁,n₂) is assumed to vary smoothly.

The ghost correction algorithm of the present invention takes advantageof the approximate shape of the phase difference function, which asdescribed above, can be obtained experimentally for test objects.Experiments indicate that variation of the function Δ (n₁,n₂) isconsiderably smaller along n₂ than n₁. In fact, for fixed n₁, variationalong n₂ is symmetric about n₂ =N/2 and Δ (n₁,n₂) can be closelyapproximated by a piecewise linear function of the form: ##EQU8##

The algorithm of the present invention takes advantage of the aboveapproximation by estimating α (n₁) and β (n₁) for each column of thedata, i.e. n₁ =0, . . . ,N_(s) -1. It consists of the following steps:

1. If raw sampled data is input, take 2-D inverse Fourier transform ofthe time data to obtain the ghosted image Y(n₁,n₂).

2. Determine signal energy for each column of the data by computing##EQU9##

3. Discard columns whose signal energy level is below a fixed threshold.This threshold is an input to the program of Appendix (A), and isdenoted by the double precision variable "snr".

4. Let S denote the set of indices of the columns whose signal energylevel is larger than the threshold snr. Estimate α (n₁) and β (n₁) forall n₁ εS. The estimation procedure will be discussed at length later.

5. Use α (n₁) and β (n₁) in equation (12) to find A(n₁,n₂) and B(n₁,n₂)for 0≦n₂ <N. From equations (9) and (10), the true object densitydistribution at (n₁,n₂) and (n₁,n₂ +N/2) are found by taking magnitudeof A(n₁,n₂) and B(n₁,n₂), respectively.

The algorithm of the present invention has been implemented in "C"programming language, and its listing is included in Appendix (A). Thename of the "C" program is "correct.c". The usage of the program isdescribed by simply invoking its executable without any parameters:

usage: correct

-g (if set, will do inverse two-dimensional FFT)

-s (snr: avoid processing columns with small signal component)

-t (ratio between--Y(n1,n2)--and--Y(n1,n2+N/2)--: must be larger forcentral columns)

-r (even to odd ratio: must be close to 1 for ghosting pixels)

-m (mse² tolerance: discard pixels whose llsq³ mse is large)

-h (ghost larger than image: 1 left, 2 right, 3 both)

-d (parameter controlling smoothness of alpha corrections)

-f (parameter controlling smoothness of beta corrections)

-a (scale to see the ghosts clearly)

-e (expand by an integer factor in each direction)

-c (convert to hexadecimal for halftone printout)

inputfile

outputfile

As it is seen, the program takes an input file representing the sampledtime data or its two-dimensional (2-D) inverse Fourier transform, andgenerates an output file containing either the binary or hexadecimalversion of the ghost-free image. Clearly, the parameters used by theprogram are set by the operator in such a way that the performance ofthe algorithm is optimized. Some of these parameters such as "-g", "-c",and "-f" are only used as switches to set flags, while others are usedto set internal variables (either integer or floating point) to specificvalues. An example of the usage of the program is as follows:

    correct -g -t100. -r1.5 -s5. -m2. -f -d.7 -b.04 -a750 -e4 inputfile outputfile

In the above example, the internal variable associated with "-t" is 100,the one associated with "-r" is 1.5, etc.

The function of various parameters used in step 1. through step 3. ofthe algorithm is as follows:

The "-g" flag determines whether or not the input ASCII file is the rawtime data or its 2-D inverse Fourier transform. Specifically, if the"-g" option is set, then the program assumes the input file to containraw time data, and computes the 2-D inverse Fourer transform of the databy successive application of the "ifft05" subroutine. This subroutinetakes inverse FFT of one-dimensional sequences, and its listing isincluded in Appendix (B).

The "-s" option sets the internal double precision variable "snr" usedin steps 2. and 3. above. A large number of the columns in the ghostedinput image correspond to empty space in the magnet, and therefore haveno signal component. To prevent the ghost correction algorithm fromprocessing these columns, the signal energy for each column of theghosted image is computed and compare it to a fixed threshold. Thisfixed threshold is the "snr" variable which is set by the "-a" option.The signal energy for the n₁ th column is defined to be: ##EQU10## whereY(n₁,n₂) denotes the value of the ghosted image at location (n₁,n₂).Thus, the AGC algorithm only processes columns whose energies exceed thethreshold set by the variable "snr". An appropriate value of "snr" is 5for input files containing raw time data, and 2×10⁸ for the onescontaining the 2-D inverse Fourier transform of the raw time data.

Step 4. of the algorithm is now described in more detail. As set forthin the Background section above, the phase difference function Δ (n₁,n₂)can be found experimentally by placing a test object in the lower orupper half of the FOV. Specifically, the ghost of a test object atlocation (n₁,n₂) with ##EQU11## appears at (n₁,n₂ +N/2), and vice versa.Thus, the location of the ghost of a pixel at (n₁,n₂) is (n₁,(n₂ +N/2)mod N). In general, if an object only fills out half of the FOV, itsghost does not overlap with itself in the reconstructed image. Underthese conditions, the phase difference function associated with theobject and the particular experimental set up can be determinedempirically for regions of the reconstructed image which correspond tothe object rather than its ghost. (Two special cases of this werediscussed in the Background section. These cases correspond to theobject being in the lower or upper half of the FOV). On the other hand,when an object fills out more than half of the FOV, its phase differencefunction can only be determined for pixels whose ghosts are notsuperimposed on pixels corresponding to other parts of the object. Thisinformation about Δ (n₁,n₂) can be exploited to find the parameters αand β as defined by equation (12). Pixels which correspond to bright(high intensity) parts whose locations correspond to either empty spacein the FOV or parts of the object with very little or no energy, will bereferred to as "ghosting pixels". Specifically, if the pixel at location(n₁,n₂) is a ghosting one, then by definition:

1. It corresponds to a high energy point in the object. Therefore

    |x(n.sub.1,n.sub.2)|>>0

2. The pixel at location (n₁,(n₂ +N/2) mod N) corresponds to a lowenergy part of an object or empty space in the FOV. That is

    |x(n.sub.1,(n.sub.2 +N/2) mod N)|≈0

Step 4 of the automatic ghost correction algorithm derives theparameters associated with the n₁ th column in two steps. Specifically,it first finds the value of the phase difference function for all the"ghosting pixels" of the column, and then solves an overdeterminedlinear system of equations to find linear least square estimates of α(n₁) and β (n₁). At this point, the key question which remains to beanswered is the way ghosting pixels are detected. The present inventionuses two criteria for classifying pixels as ghosting ones.

The first criterion is a direct consequence of equation (11) and thesecond part of the definition of ghosting pixels. It takes advantage ofthe fact that the magnitude of the even and parts of a ghosting pixel atlocation (n₁,n₂) are identical. Specifically, if Y(n₁,n₂) is a ghostingpixel, from equation (11):

    |Y.sub.even (n.sub.1,n.sub.2)|=|Y.sub.odd (n.sub.1,n.sub.2)|=|A(n.sub.1,n.sub.2)|=x(n.sub.1,n.sub.2)                                                (14)

Thus, if the ratio between the magnitude of even and odd parts of thepixel (the "eoratio") at location (n₁,n₂) are equal or somewhat close toeach other, then Y(n₁,n₂) can be classified as a ghosting pixel. In theprogram listing of the automatic ghost correction algorithm in Appendix(A), if the ratio between the magnitudes of even and odd parts of apixel are between [1/eoratio, eoratio], then the pixel is classified asa ghosting one. As seen, "eoratio" denotes a double precision variablewhich is input to the program by the operator via parameter -r. It isimportant to note that under this criterion, the conditions for Y(n₁,n₂)and Y(n₁,n₂ +N/2) to qualify as ghosting pixels are identical Thus, ifequation (14) is exactly or approximately satisfied, then there is anambiguity as to which pixel is the ghosting one. As discussed below, thesecond criteria for ghosting pixel detection helps to resolve thisambiguity.

The appropriate values of "eoratio" is anywhere between 1 and 2. In mostof the examples of the next section, the "-r" option (the internalvariable "eoratio") is set to 1.5. If it is set to small values (i.e.too close to 1), then the number of ghosting pixels will be too small,and therefore estimation of α (n₁) and β (n₁) will not be robust. On theother hand, if it is set to larger values such as 2 or even 3, then thechosen pixels might not necessarily be the ghosting ones. This willincrease the error in the observations for the linear least-squaresestimation, and therefore will result in less accurate estimation of α(n₁) and β (n₁).

The second criteria takes advantage of the definition of ghostingpixels. To describe this condition, equations (1) and (2) are rewrittenin the following way: ##EQU12## As expected, if θ (n₁,n₂)=φ (n₁,n₂)

or equivalently

Δ (n₁,n₂)=0

the observed image, Y(n₁,n₂), becomes ghost free and is identical to thetrue object density function x(n₁,n₂). Recall that if the pixel atlocation (n₁,n₂) is a ghosting one, then by definition, the magnitude ofx(n₁,n₂) must be large (i.e. not at the noise level) and the magnitudeof x(n₁,n₂ +N/2) must be very small (i.e. not at the noise level). Fromequations (15) and (16), if the difference between θ (n₁,n₂) and φ(n₁,n₂) is small (or equivalently Δ (n₁,n₂) is small), then a ghostingpixel at (n₁,n₂) results in a large value of the following ratio:##EQU13##

Specifically, as Δ (n₁,n₂) changes from 0 to π, the above ratio changesfrom ∞ to 0. It has been found experimentally that Δ (n₁,n₂) is smallerfor columns closer to the center of the magnet (i.e. around n₁ =N_(s)/2). This implies that the ratio of equation (17) becomes larger as theindex of the column under investigation becomes closer to N_(s) /2.Thus, the second criterion for detecting ghosting pixels of the n₁ thcolumn consists of

1. Computing the quantity shown in equation (17) for each pixel.

2. Comparing this ratio to a fixed threshold associated with thatcolumn.

Clearly, this threshold is column dependent and becomes larger as theindices of the column become closer to N_(s) /2. In the program listedin Appendix (A), the threshold for the n₁ th column is:

the internal variable "threshold" set by "-t" option, if n₁ =N_(s) /2.

decreases linearly with |n₁ -N_(s) /2|for |n₁ -N_(s) /2|<15.

is equal to 1 for |n₁ -N_(s) /2|>15.

In most of the following examples, the "-t" option (the internal doubleprecision "threshold") is set to 100. In general, the appropriate valuefor the "-t" option depends on the amount of ghosting in the centralpart of the image, and lies between 1 and 10000. If the reconstructedimage suffers from considerable amount of ghosting in the centralcolumns, then "threshold" must be set at a small value e.g. 1.Otherwise, it should be set at a larger value, say 100, so that thecenter 30 columns of the image remain more or less unchanged by thealgorithm. If the central columns of an image are ghost-free, settingthe "-t" option at small values might result in unnecessary distortionsin these columns.

To summarize, the first ghost detection criterion checks the ratiobetween the magnitudes of even and odd parts of the pixel at location(n₁,n₂). If this ratio is close to one, then either Y(n₁,n₂) or Y(n₁,n₂+N/2) are classified as ghosting pixels. To resolve this ambiguity andto improve the detection procedure, a second criterion is used whichcomputes the ratio shown in equation (17). For columns close to thecenter of the magnet, large values of this ratio imply a ghosting pixelat location (n₁,n₂). However, for columns further away from the center,the second criterion becomes more or less inconclusive, and other waysmust be found to overcome the ambiguity problem of the first criterion.The present invention uses the apriori knowledge about the approximateshape of the phase difference function in order to resolve thisambiguity. Detailed experimental procedures for obtaining the phasedifference function was described in the Background section above.Unlike that "look up" table approach, the present invention does notrequire detailed and exact values of the phase difference function. Infact, the algorithm of the present invention only needs to know as muchabout Δ (n₁,n₂) as to make binary decisions.

From the description of the automatic ghost correction algorithm, it isclear that the ghosting pixel detection part is rather hueristic. Todecrease the sensitivity of the algorithm to this part, and to improvethe estimation of the phase difference function, the following measuresare preferably employed:

From classical results in estimation theory, it is clear that the errorin estimating the parameters of Δ (n₁,n₂), i.e. α and β, is reduced asthe number of observations is increased. In this case, the observationsare the empirical value of the phase difference function for theghosting pixels. Experimental results seem to indicate that for columnswhose number of observations (or equivalently the number of ghostingpixels) is small, the error in estimating β in equation (12) becomesrather large. This error manifests itself as large magnitude for β,resulting in unrealistic values of the phase difference function. Sinceexperimental results indicate that the magnitude of β is small for mostcolumns, the present invention sets to β=0 for columns whose number ofghosting pixels is less than a fixed integer. In the program listing ofAppendix (A), this integer has been chosen to be 8. To make theestimation part of the algorithm even more robust, β is set to zero whenthe magnitude of its estimated value exceeds a certain threshold. Theoptimal value for this threshold was found empirically from theapproximate shape of the phase difference function for various testobjects. For the program listed in Appendix (A), this threshold was setto 0.08.

A second measure taken to improve the robustness of the algorithm is todiscard ghosting pixels whose least-square residue is too large.Specifically, if i₁,i₂, . . . , i_(upper) <N/2and j₁,j₂, . . . ,j_(lower) ≧N/2 denote the indices of the ghosting pixels of the n₁ thcolumn, taking into account equation (12), the linear least-squaresestimation of α and β consists of solving the following overdeterminedlinear system of equations: ##EQU14##

If α (n₁) and β (n₁) denote the solution of the above linearleast-squares problem, the residue of the ghosting pixels at (n₁,i₁) and(n₁,j₁) are defined to be: ##EQU15## and the means square error isdefined to be ##EQU16##

To reduce the likelihood of false alarm (i.e. declaring non-ghostingpixels as ghosting ones), pixels for which the ratio between theirresidue and the mean squared error exceeds a predetermined. threshold"mse" are discarded.

The "-m" option sets the threshold "mse". Appropriate values for the"-m" option can range anywhere between 1 and 3. Clearly, if "mse" is toosmall, then most of the already detected ghosting pixels will bediscarded for the second estimation process. On the other hand, if "mse"is too large, picels which have been misclassified as ghosting ones arenot discarded and therefore result in large amounts of error inobservations used for the linear least-squares estimation process. Inmost of the examples of the next section, the "-m" option is set to 2.

Another measure to improve the robustness of the algorithm is the -hoption. Recall that the ghost detection part of the algorithm firstchecks the ratio between the magnitudes of even and odd parts of thepixel at location (n₁,n₂). If this ration is close to 1 (or morespecifically is in the range [1/eoratio, eoratio]) then either Y(n₁,n₂)or Y(n₁,(n₂ +N/2) mod N) are classified as ghosting pixels. If the phasedifference function is less than π then the magnitude of the ghost (i.e.the ghosted pixel) is smaller than that of the object (i.e. the ghostingpixel). On the other hand, if the phase difference function is largerthan π, the magnitude of the ghost becomes larger than that of theobject causing it. Since we expect the phase difference function to besmall (less than π) in the center of the magnet, the differentiationbetween ghosting and ghosted pixels is straightforward for the centralcolumns of the image, and therefore there is not any ambiguity incomputing the phase difference function. In fact, as we mentionedearlier, the "-t" option takes advantage of this in order todetect/differentiate ghosting and ghosted pixels for the central 30columns. Specifically, if the ratio shown in equation (17) is largerthan the "eoratio" parameter set by the "-t" option, then Y(n₁,n₂) isthe ghosting pixel and Y(n₁,n₂ +N/2) is the ghosted pixel. On the otherhand if the ratio of equation (17) is smaller than 1/eoratio, thenY(n₁,n₂ +N/2) is the ghosting and Y(n₁,n₂) is the ghosted pixel.However, this clear distinction between ghosting and ghosted pixels ispossible only if the phase difference function is known to be less thanπ, which in out situation corresponds to the central columns of theimage. For columns away from the center, the phase difference functionmight become larger than π, thus creating an ambiguity about theghosting and ghosted pixels. To resolve this ambiguity, the operator hasto provide the algorithm with a clue as to whether or not the phasedifference function is larger than π. Fortunately, this is a simplevisual task since in the areas of image with large phase differencefunction (i.e. larger than π) the ghosts are brighter than the actualobject causing the ghost. Therefore, if the ghost has larger intensitythan the object in the areas to the left of central columns, "-h" optionmust be set to 1. Similarly, if the ghost has larger intensity than theobject in the areas to the right of central columns, "-h" option must beset to 3. As shown below in most imaging situations, the ghost has lowerintensity than the actual object itself, and thus there is no need toset the "-h" option to any value. However, in few of the heart imagesshown in FIGS. 4-11, the ghosts become larger than the objects and the"-h" option is needed to guide the program to correct them.

The "-d" option sets the internal double precision parameter "smooth"which smoothes the values of for neighboring columns. Specifically, ifthe value of for the n₁ th column is different from those of the pastfour columns by more than an amount specified by "smooth", then α (n₁)and β (n₁) are discarded and the n₁,th column appears unchanged in theoutput image Appropriate values for "smooth" can range between 0.3 and2. Clearly, if "-d" option is set to a small value, e.g. 0.3, then α andβ for most columns will be discarded and the final output will resemblethe ghosted input image to a great extent. On the other hand, if "-d"option is set to a large value like 2, then the amount of smoothing isminimized, and the final image might have some columns which stand outamont their neighboring columns because of their drastically, differentvalues of α and β.

The "-b" option sets the internal double precision parameter "betasmth"which smoothes the values of β for neighboring columns. The functionaldescription of this variable is similar to that of the "-d" option.Appropriate values for "betasmth" however, can range between 0.01 and 1.

The "-f" flag determines whether or not the processed ghost-free imageneeds to be rearranged so that the center of the magnet coincides withthe center of the image. This flag is set for all the examples shown inthe next section.

The "-a" option sets the internal double variable "scale" which scalesthe processed image before it is written in the output file. Thisparameter is normally set anywhere between 500 and 1000.

The "-e" option sets the internal variable "expand" which expands thefinal output image by an integer in each direction. The expansion isbasically done by repeating each pixel a fixed number of times in x andy directions.

The "-c" flag determines whether the processed image is written inbinary or hexadecimal format. The hexadecimal format is necessary forgenerating halftone images with PostScript commands and the AppleLaserWriter.

"Inputfile" denotes the name of the input file which contains the ASCIIcharacters representing the time raw data or its two-dimensional inverseFourier transform. If such an input file does not exist, the programexists while printing "input file does not exist".

"Outputfile" denotes the name of the output file which contains numbersbetween 0 and 255 representing the processed image. The format of thesenumbers is either in binary or hexadecimal depending on whether or notthe "-c" flag is set.

Examples of images processed by the ghost elimination algorithm areshown in FIGS. 4-10. For each example, the ghosted image, the processedghost-free image, and the parameters used with the algorithm will beshown.

FIG. 4a shows a ghosted heart image which was processed by the ghostelimination algorithm with the following parameters:

    correct -g -t100. -r1.5 -s5. -m3. -h1 -f -d2. -b1. -a1000 -e4 htnum htproc

The processed ghost-free version is shown in FIG. 4b. Clearly, the AGCalgorithm has done an excellent job of removing the ghosts. As it isseen, the magnitude of the ghost in the left side of FIG. 4a is largerthan that of the object causing it. Therefore, the parameter "-h" had tobe set at 1 in order to remove the ambiguity in the phase differencefunction for the columns to the left of the image.

The second example of the ghost elimination algorithm is shown in FIG.5. The original ghosted heart image is shown in FIG. 5a, and itsprocessed ghost-free version is shown in FIG. 5b. The parameters of thealgorithm were:

    correct -g -t100. -r1.5 -s5. -m2. -f -d1. -b.05 -a1000 -e4 heart5num heart5proc

A third example of the algorithm is shown in FIG. 6a, and its processedghost-free version is shown in FIG. 6b. The parameters of the algorithmwere:

    correct -g -t100. -r1.5 -s5. -m2. -f -d2. -b.1 -h1 -a750 -e4 heart4num heart4proc

Similar to FIG. 4a, the magnitude of the ghost in the left side of FIG.6a is larger than that of the object causing it. Thus, "-h" option hadto be set to 1, to remove the ambiguity of the phase difference functionfor the columns to the left of the ghosted image.

A fourth example of the algorithm is shown in FIG. 7. The originalghosted heart image is shown in FIG. 7a, and its processed ghost-freeversion is shown in FIG. 7b. The parameters of the algorithm were:

    correct -g -t100. -r1.5 -s5. -m2. -f -d.7 -b.03 -h3 -a500 -e4 heart3num heart3proc

As it is seen in FIG. 7a, the magnitude of the ghost in the right andleft part of the ghosted image is slightly larger than that of theobject causing it. Therefore, the "-h" option is set at 3, to remove theambiguity associated with the phase difference function.

A fifth example of the algorithm is shown in FIG. 8. The originalghosted heart image is shown in FIG. 8a, and its processed ghost-freeversion is shown in FIG. 8b. The parameters of the algorithm were:

    correct -t10. -r1.3 -s200000000. -m2. -f -d.5 -b1. -a1000 -e4 newhtnum newhtproc

Note that in the above example, the input file contained the 2-D inverseFourier transform of the raw time date. Therefore, the "-s" option hadto be set at a different value from the previous examples.

A sixth example of the algorithm is shown in FIG. 9. The originalghosted heart image is shown in FIG. 9a, and its processed ghost-freeversion is shown in FIG. 9b. The parameters of the algorithm were:

    correct -g -t1. -r1.5 -s5. -m3. -f -d.6 -b.06 -a1000 e4 lvrnum lvrproc

A seventh example of the algorithm of the algorithm is shown in FIG. 10.The original ghosted heart image is shown in FIG. 10a, and its processedghost-free version is shown in FIG. 10b. The parameters of the algorithmwere:

    correct -g -t100. -r1.5 -s5. -m2. -f -d1. -b.1 -a1000 -e4 body2num body2proc

An eighth example of the algorithm is shown in FIG. 11. The originalghosted heart image is shown in FIG. 11a, and its processed ghost-freeversion is shown in FIG. 11b. The parameters of the algorithm were:

    correct -g -t100. -r1.5 -s5. -m3. -f -d.7 -b.06 -a1000 -e4 body3num body3proc

Although the present invention has been described in connection with apreferred embodiment thereof, many other variations and modificationswill now become apparent to those skilled in the art without departingform the scope of the invention. It is preferred, therefore, that thepresent invention be limited not by the specific disclosure herein, butonly by the appended claims ##SPC1##

What is claimed is:
 1. A method of cancelling .[.ghosts.]. .Iadd.a ghost.Iaddend.from .Iadd.a .Iaddend.NMR .[.images.]. .Iadd.image.Iaddend.,comprising the steps of:(a) taking a two dimensional inverse Fouriertransform of a raw NMR signal to obtain a ghosted image Y(n₁,n₂); (b)computing the signal energy for each column of said Fourier transformedsignal using ##EQU17## (c) discarding the columns whose signal energylevel are below a predetermined threshold; (d) .[.estimating α (n₁) andβ (n₁) for.]. .Iadd.identifying ghosting pixels in .Iaddend.eachremaining column of data; i.e. n₁ =0, . . . N_(s) -1, .[.by:.]..[.(i).]. .Iadd.(e) .Iaddend.finding the .Iadd.actual .Iaddend.phasedifference .[ƒ]. Δ.Iadd._(actual) .Iaddend.(n₁,n₂) for .[.all.]..Iadd.the .Iaddend.ghosting pixels .[.of the column; and.]. .[.(ii).]..Iadd.(f) .Iaddend.solving the following simultaneous equations to findlinear least square estimates of α (n₁) and β (n₁): ##EQU18## .[.(e).]..Iadd.(g) .Iaddend.using .Iadd.the linear least square estimates of.Iaddend.α (n₁) and β (n₁) .[.in the above equation.]. to find .[.the.]..Iadd.an estimated .Iaddend.phase difference .Iadd.function .Iaddend.Δ(n₁,n₂) for 0≦n₂ <N.[.; and.]. .Iadd.using .Iaddend. ##EQU19## .[.(f).]..Iadd.(h) .Iaddend.using .Iadd.the estimated phase difference function.Iaddend.Δ (n₁,n₂) in ##EQU20## to find A(n₁,n₂) and B(n₁,n₂) for 0≦n₂<N, where the dimensions of the .[.reconstructed.]. image .Iadd.datax(n₁,n₂) corresponding to the NMR image .Iaddend.are N×N_(s) .Iadd.andwhere .Iaddend. ##EQU21##
 2. The method of claim 1, wherein the ratio ofthe magnitude of even and odd parts of a pixel at location (n₁,n₂) isused to determine whether a pixel is a true image or a ghost.
 3. Themethod of claim 2, wherein the pixel is classified as a ghost if saidratio is approximately equal to one.
 4. The method of claim 2, whereinthe following ratio is computed, indicating that a pixel near the centerof the magnetic cord is a ghost when the ratio extends to apredetermined value.
 5. The method of claim 1, wherein columns withsmall signal components are not processed.
 6. The method of claim 1,wherein pixels at location (n₁, n₂) whose linear base squares meansquare error is larger than a predetermined amount are discarded whendetermining the phase difference function. .Iadd.7. The method of claim1 wherein the actual phase difference Δ_(actual) (n₁,n₂) is calculatedusing:Δ_(actual) (n₁,n₂)=phase (Y_(even) (n₁,n₂))-phase (Y_(odd)(n₁,n₂))..Iaddend..Iadd.8. A method of producing a magnetic resonanceimage, comprising the steps of:applying a bipolar readout gradientfield; sampling NMR data by scanning N lines of the k_(x) -k_(y) spacein connection with the application of the bipolar readout gradientfield, with reversal of the direction of data sampling on odd and evenlines; processing the NMR data to generate image data representative ofa magnetic resonance image such that there is substantially no phasedifference between even and odd parts of the imagedata..Iaddend..Iadd.9. The method of claim 1 wherein the phasedifference between even and odd parts of the transform data is afunction that varies in at least one of the x and ydirections..Iaddend..Iadd.10. The method of claim 8 wherein the step ofsampling NMR data comprises introducing a time delay in sampling thedata, the time delay being selected to provide substantially no phasedifference between even and odd parts of the imagedata..Iaddend..Iadd.11. The method of claim 10 wherein the step ofsampling NMR data comprises introducing a time delay in sampling thedata, the time delay being selected to provide substantially no firstorder phase difference between even and odd parts of the imagedata..Iaddend..Iadd.12. The method of claim 8 wherein the step ofprocessing the NMR data comprises transforming the NMR data intotransform data..Iaddend..Iadd.13. The method of claim 12 wherein thestep of processing the NMR data comprises determining a functionrepresentative of the phase difference between even and odd parts of thetransform data..Iaddend..Iadd.14. The method of claim 12 wherein thestep of transforming the NMR data comprises taking an inverse Fouriertransform of the NMR data..Iaddend..Iadd.15. The method of claim 14wherein the step of transforming the NMR data comprises taking atwo-dimensional inverse Fourier transform of the NMRdata..Iaddend..Iadd.16. The method of claim 12 wherein the transformdata are denoted by Y(n₁,n₂) and where: ##EQU22##.Iadd.17. The method ofclaim 16 further comprising the step of computing the signal energy foreach column of the transform data using .Iadd.18. The method of claim 17further comprising the step of discarding columns with a signal energybelow a predetermined threshold..Iaddend..Iadd.19. The method of claim16 wherein there is an actual phase difference Δ_(actual) (n₁,n₂)between Y_(even) (n₁,n₂) and Y_(odd) (n₁,n₂)..Iaddend..Iadd.20. Themethod of claim 19 wherein the actual phase difference Δ_(actual)(n₁,n₂) is calculated using: Δ_(actual) (n₁,n₂)=phase (Y_(even)(n₁,n₂))-phase (Y_(odd) (n₁,n₂))..Iaddend..Iadd.21. The method of claim19 wherein the step of processing the NMR data comprises using thetransform data Y(n₁,n₂) and the actual phase difference Δ_(actual)(n₁,n₂) to calculate the image data, wherein the image data are denotedby x(n₁,n₂)..Iaddend..Iadd.22. The method of claim 21 wherein the stepof processing the NMR data comprises estimating a phase differencefunction Δ(n₁,n₂) from the actual phase difference Δ_(actual)(n₁,n₂)..Iaddend..Iadd.23. The method of claim 22 wherein the estimatedphase difference function Δ(n₁,n₂) is determined by(i) calculating aplurality of values of the actual phase difference Δ_(actual) (n₁,n₂);(ii) using the plurality of values of the actual phase differenceΔ_(actual) (n₁,n₂) to find linear least square estimates of α(n₁) andβ(n₁) using ##EQU23## (ii) using the linear least square estimates ofα(n₁) and β(n₁) to find the estimated phase difference function Δ(n₁,n₂)for 0≦n₂ <N using ##EQU24## .Iadd.24. The method of claim 22 wherein thestep of processing the NMR data comprises using the transform dataY(n₁,n₂) and the estimated phase difference function Δ(n₁,n₂) tocalculate the image data x(n₁,n₂)..Iaddend..Iadd.25. The method of claim24 wherein the image data x(n₁,n₂) are calculated by: (i) using theestimated phase difference function Δ(n₁,n₂), Y_(even) (n₁,n₂), andY_(odd) (n₁,n₂) to find A(n₁,n₂) and B(n₁,n₂) using ##EQU25## where##EQU26## (ii) using the respective magnitude of A(n₁,n₂) and B(n₁,n₂)to find x(n₁,n₂) and x(n₁,n₂ +N/2)..Iaddend..Iadd.26. A method ofreducing a ghost from a NMR image, comprising the steps of:(a) taking atwo dimensional inverse Fourier transform of a raw NMR signal to obtaina ghosted image Y(n₁,n₂); (b) identifying ghosting pixels in columns ofY(n₁,n₂); i.e. n₁ =0 . . . N_(s) -1, (c) finding the actual phasedifference Δ_(actual) (n₁,n₂) for the ghosting pixels (d) solving thefollowing simultaneous equations to find linear least square estimatesof α(n₁) and β(n₁); ##EQU27## (e) using the linear least squareestimates of α(n₁) and β(n₁) to find an estimated phase differencefunction Δ(n₁,n₂) for 0≦n₂ <N using ##EQU28## (f) using the estimatedphase difference function Δ(n₁,n₂) in ##EQU29## to find A(n₁,n₂) andB(n₁,n₂) for 0≦n₂ <N, where the dimensions of the image data x(n₁,n₂)corresponding to the NMR image are N×N_(s) and where ##EQU30## .Iadd.27.A method of reducing a ghost from a NMR image, comprising the steps of:(a) taking a two dimensional inverse Fourier transform of a raw NMRsignal to obtain a ghosted image Y(n₁,n₂);(b) calculating Y_(even)(n₁,n₂) and Y_(odd) (n₁,n₂) using ##EQU31## (c) estimating a phasedifference function Δ(n₁,n₂) between Y_(even) (n₁,n₂) and Y_(odd)(n₁,n₂); and (d) using the estimated phase difference function Δ(n₁,n₂)in ##EQU32## to find A(n₁,n₂) and B(n₁,n₂) for 0≦n₂ <N, where the imagedata corresponding to the NMR image are denoted by x(n₁,n₂), and where##EQU33## .Iadd.28. The method of claim 27 wherein the estimated phasedifference function is predetermined..Iaddend..Iadd.29. The method ofclaim 28 wherein the estimated phase difference function is estimatedfrom at least a portion of a Fourier transform of the NMRsignal..Iaddend..Iadd.30. The method of claim 29 wherein the estimatedphase difference function is estimated from the ghosted imageY(n₁,n₂)..Iaddend.